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ABSTRACT 

MIT  Lincoln  Laboratory  (MIT/LL)  has  been  developing  a  computational  method  to  predict  the 
maximum  structural  response  due  to  a  pyrotechnic  shock  input  using  Statistical  Energy  Analysis 
(SEA).  It  had  been  previously  understood  that  since  the  pyrotechnic  shock  is  not  a  steady  state 
event,  traditional  SEA  method  may  not  applicable.  A  new  analysis  methodology  effectively 
utilizes  the  traditional  SEA  modeling  by  employing  an  energy  constraint  using  the  Shock 
Response  Spectrum  (SRS)  as  the  steady  state  input  condition.  Responses  are  recovered  in  the 
SRS  and  can  be  used  to  develop  pyrotechnic  shock  design  and  test  requirements  for  various 
subsystems  of  the  model.  Response  accelerations  in  the  time  domain  corresponding  to  the 
response  spectra  are  estimated  using  the  Peak  Ratio  (PR).  The  response  waveform  is  obtained 
using  the  inverse  transform  of  the  product  of  the  input  acceleration  and  the  transfer  function  of 
the  SEA  model.  The  most  probable  input  acceleration  may  be  obtained  using  the  stochastic  SRS 
decomposition  method.  Response  waveforms  can  be  further  refined  based  on  phase  perturbations 
using  the  wave  numbers  of  the  SEA  subsystems.  A  couple  of  validation  examples  are  provided  to 
demonstrate  the  new  approach. 

KEY  WORDS:  Peak  Ratio,  phase  perturbation,  pyrotechnic  shock,  shock  response  spectrum. 
Statistical  Energy  Analysis. 

SYMBOLS 

g  =  gravitational  constant,  9.81  m/s^ 

gpk  =  temporal  peak  acceleration  in  gravitational  constant 

g^  =  root  mean  square  acceleration  in  gravitational  constant 

INTRODUCTION 

Launch  vehicles  and  space  systems  use  various  types  of  pyrotechnic  shock  devices  during  their 
mission  lives.  While  significant  progress  has  been  made  in  the  last  few  decades,  high  fidelity 
simulation  of  structural  responses  due  to  a  pyrotechnic  shock  input  is  still  one  of  the  challenging 
problems  in  space  technologies.  As  noted  by  NASA^  and  ESA^  references,  the  current  shock 
prediction  methods  do  not  consistently  produce  successful  results.  Analysts  frequently  utilize  the 
Einite  Element  Method  (FEM)  using  the  modal  transient  or  direct  integration  methods  to  predict 
structural  responses  due  to  a  pyrotechnic  shock  input.  Unfortunately,  there  are  many  problems 
associated  with  this  approach.  First,  pyrotechnic  shock  environments  are  broadband  and  typically 
analyzed  at  least  to  10000  Hz.  At  low  frequency,  structural  modes  may  be  predicted  using  a 
deterministic  method  such  as  the  FEM,  however  at  high  frequency  structures  typically  produce 
modally  dense  responses  with  significant  overlaps.  Since  the  responses  at  high  frequency 
typically  converge  to  the  mean  value,  an  alternative  method  that  can  simulate  an  ensemble 
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average  response  would  be  better  suited.  Second,  in  order  to  use  the  FEM  at  high  frequency,  the 
model  would  require  a  large  number  of  elements  because  it  is  typical  to  use  at  least  4  to  6 
elements  per  each  wavelength.  Therefore,  it  becomes  impractical  to  use  the  FEM  for  a  full 
system  analysis  at  high  frequency.  In  addition  to  the  large  model  size,  high  frequency  responses 
would  be  severely  influenced  by  the  complex  modal  interactions  among  various  structural 
members  such  as  intervening  structures,  interfaces,  junctions,  etc.  As  a  result,  it  becomes  too 
complex  to  track  the  response  deterministically  at  high  frequency.  Therefore,  FEM  modeling  is 
typically  applicable  to  a  simple  structure  when  the  input  is  a  pyrotechnic  shock  related  event. 

In  this  paper.  Statistical  Energy  Analysis  (SEA)  is  used  as  an  alternative  method  to  produce  a 
response  at  high  frequency.  Since  the  SEA  was  originally  developed  to  produce  an  ensemble 
average  response  at  high  frequency^,  it  would  be  better  suited  for  pyrotechnic  shock  response 
predictions.  However,  most  of  the  times,  the  method  is  simplified  for  steady  state  problems  and 
as  a  result,  traditional  SEA  modeling  may  not  be  directly  applicable  to  pyrotechnic  shock 
problems.  Manning"^  investigated  a  structural  response  of  a  simple  model  subjected  to  a  transient 
input  forcing  function  using  the  SEA  modeling.  In  his  approach,  power  balance  equation  was  not 
simplified  for  steady  state  cases  and  the  equation  was  solved  at  each  time  step  with  the  steady 
state  assumption  for  the  coupling  loss  factors.  This  approach  is  frequently  referenced  as  the 
transient  SEA  method  in  the  literature.  Dalton^  used  the  FEM  and  SEA  modeling  to  produce  a 
response  for  a  broad  frequency  band  of  interest.  Input  forcing  function  was  estimated  using  the 
CTH  Hydrocode  developed  by  the  Sandia  National  Laboratories^.  This  approach  is  often  cited 
Virtual  Mode  Synthesis  and  Simulation  (VMSS)^  because  it  is  based  on  curve  fitting  of  the 
maximum  magnitudes  of  imaginary  single  degrees  of  freedom  system  models  to  the  Frequency 
Response  Function  (FRF)  of  the  SEA  model.  Based  on  improved  mid  frequency  predictions^, 
Borello^  developed  a  shock  prediction  method  using  the  SEA  modeling.  In  his  approach,  modal 
data  (either  from  a  FEM  model  or  a  test)  are  used  with  the  SEA  model  to  produce  a  response 
waveform  with  the  random  phase  distribution.  Input  forcing  functions  simulating  several  types  of 
pyrotechnic  shock  sources  have  been  developed  in  the  wave  number  domain  as  a  result. 

In  the  upcoming  sections,  a  new  pyrotechnic  shock  analysis  methodology  using  the  SEA 
modeling  is  introduced.  A  numerical  example  and  a  preliminary  study  of  the  complex  validation 
model  are  provided  to  demonstrate  the  new  approach. 

PROPOSED  PYROTECHNIC  SHOCK  ANALYSIS  METHODOLOGY  USING  SEA 

Figure  1  shows  the  overview  of  the  proposed  method  using  the  SEA  modeling.  The  new 
approach  differs  from  the  other  methods  in  several  categories.  First,  unlike  the  other  methods, 
the  proposed  method  does  not  require  input  forcing  functions.  Instead,  the  new  approach  uses  the 
SRS  as  an  energy  constraint  condition  to  the  SEA  model.  The  SRS  is  the  maximum  absolute 
temporal  response  of  a  single  degree  of  freedom  system  due  to  an  arbitrary  base  driven  transient 
input  at  each  prescribed  frequency.  For  a  given  transient  waveform,  only  one  SRS  is  available. 
Alternatively,  the  SRS  can  be  interpreted  as  the  maximum  response  state  of  a  system  that  does 
not  change  in  time.  This  attribute  allows  the  use  of  the  traditional  SEA  modeling.  By  using  the 
SRS  as  the  input  condition,  coupling  loss  factors  of  the  SEA  model  would  simulate  energy 
transfers  among  subsystems  at  the  maximum  response  state.  Since  the  pyrotechnic  shock  design 
and  test  requirements  are  provided  in  the  SRS,  this  approach  can  be  easily  implemented  as  the 
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energy  constraint  condition  at  the  interface  of  the  SEA  model  to  develop  pyrotechnic  shock 
design  and  test  requirements  for  various  subsystems  of  the  model. 


Figure  1:  Overview  of  the  proposed  shock  analysis  method  using  SEA. 

Second,  since  the  SRS  is  used  as  the  input,  responses  are  recovered  in  the  SRS  as  well.  Estimated 
response  acceleration  in  the  time  domain  at  each  frequency  can  be  developed  using  the  Peak 
Ratio  (PR)*°.  The  PR  describes  the  relationship  between  the  input  in  the  time  domain  and  the 
response  spectra  in  the  frequency  domain.  In  reality,  the  relationship  between  the  input  and  the 
response  can  be  described  using  the  combination  of  the  PR,  Energy  Ratio  (ER),  and  the  phase^°. 
However,  among  the  three  variables,  the  PR  is  the  most  dominant,  and  in  this  paper,  only  the  PR 
is  used  to  estimate  the  response  acceleration  in  the  time  domain  at  each  frequency.  Equations  1 
and  2  describe  the  PR  and  the  ER,  respectively.  The  obtained  response  acceleration  magnitudes 
are  RSS’ed  to  produce  the  maximum  global  response  magnitude. 

PR=—  (1) 

Spk 

ER  =  —  (2) 

^  tins 

Third,  optionally,  a  response  waveform  may  be  estimated  using  the  Transfer  Eunction  (TE)  of  the 
SEA  model.  Most  of  the  times,  the  input  acceleration  is  not  available  and  in  this  case  the  SRS 
requirement  may  be  decomposed  into  the  most  probable  input  accelerations  using  the  stochastic 
SRS  decomposition  method”.  Inverse  Fourier  transform  of  the  product  of  the  Fourier  transform 
of  the  input  acceleration  and  the  TF  of  the  SEA  model  produces  the  response  waveform. 
Magnitude  of  the  response  waveform  would  likely  be  close  to  the  maximum  global  response 
magnitude  obtained  from  the  SRS  constraint  condition.  However,  in  case  there  is  a  difference, 
response  waveform  may  be  scaled  to  match  the  maximum  global  response  magnitude  obtained 
from  the  SRS  energy  constraint  condition.  In  this  proposed  method,  the  response  waveform 


29th  Aerospace  Testing  Seminar,  October  2015 


would  carry  the  phase  of  the  input  acceleration  to  the  response  location.  This  would  be  a  largely 
incorrect  presumption  and  a  eoneeptual  mathematical  model  that  eould  be  used  to  simulate  the 
phase  change  (between  the  source  and  the  receiver)  has  been  developed  based  on  wave  numbers 
of  the  SEA  subsystems.  Three  types  of  waves  ineluding  flexural,  longitudinal,  and  shear  could  be 
used  in  the  formulation  to  perturb  the  phase  of  the  response  waveform.  The  conceptual  model 
also  allows  a  weighting  eould  be  added  to  the  perturbation  based  on  modal  energy  or  total  energy 
of  the  SEA  model.  Figure  2  shows  the  mathematical  concept  of  the  phase  perturbation  method. 


n_subnj 

n  modfX\j>« 

=  Phase(  ■  X  ) 

}  mr  l 


w  =  frequency 

d(uj)  =  phase  shift  between  input  and  response 
locations  at  frequency  u> 

njsubsys  =  number  of  subsystems  in  SEA  model 
jjss  =  jth  subsystem  in  SEA  model 
n_modetype  =  number  of  different  types  of  wave 
modes 

j_mt  =  f'  type  of  wave  mode 
W  =  weighting  factor  used  to  represent  relative 
differences  in  acceleration  amplitude  due  to 
contributions  from  different  types  of  wave  modes 
k  =  wave  number 

C  -  distance  wave  travels  through  subsystems 


Figure  2:  Conceptual  phase  perturbation  mathematical  model. 

NUMERICAL  EXAMPLE 

Figure  3  shows  a  numerical  validation  model  comprised  of  five  beam  and  three  plate  SEA 
subsystems.  The  SRS  input  was  applied  as  the  energy  constraint  condition  normal  to  the  Small 
Side  Plate  subsystem.  Figure  3  also  shows  the  comparable  FEM  model  developed  with  48  CBAR 
elements  for  the  Beam,  23104  QUAD4  elements  for  each  Large  Plate,  24  CBAR  elements  for 
each  Small  Beam,  and  1742  QUAD4  elements  for  the  Small  Side  Plate.  Figure  4  shows  the 
transient  input  acceleration  applied  in  normal  direction  at  the  center  of  the  Small  Side  Plate  for 
the  modal  transient  analysis.  Information  of  the  validation  model  can  be  found  in  Table  1. 


Table  1;  Information  of  the  validation  model. 


Matertal 

Size 

Mass 

Boundary 

Condition 

Damping 

Beam 

Aluminum 

0.036m  X  O.D36m  X 
O.G09m 

2.121  Kg 

Free  top 

Clamped  to  plate 

Q=10 

Small  Beam  (4) 

Aluminum 

0.036m  X  O.D36m  X 
0,305m 

1.061  Kg 

Clamped 

Q=10 

Large  Plate  (2} 

Aluminum 

1,026m  X  1.926m  X 

0,010m 

95.400Kg 

Clamped 

Q=1D 

Small  Side  Plate 

Aluminum 

0.3S5m  X  0.385m  X 

0,010m 

3.816Kg 

Clamped  to  four 
beams 

Q=1D 
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Figure  3:  Numerical  validation  model  (SEA  and  FEM). 
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Figure  4:  Pyrotechnic  shock  transient  input  acceleration. 


Figure  5  shows  the  response  measured  at  Node  12359  of  the  Top  Plate  using  the  modal  transient 
method.  Temporal  peak  acceleration  values  of  3Ag  and  -3.6g  have  been  observed  from  the 
response.  Figure  6  shows  responses  of  3  SEA  subsystems  recovered  in  the  SRS  with  the  SRS 
input  applied  as  the  energy  constraint  normal  to  the  Small  Side  Plate  subsystem.  The  SRS  input 
has  been  derived  from  the  transient  input  acceleration  in  Figure  4.  The  response  spectra  in  Figure 
6  are  divided  by  the  3  PR  values  (3,  5,  7)  in  Figure  7  to  produce  the  estimated  response 
acceleration  in  the  time  domain  at  each  frequency.  The  maximum  global  response  magnitude  is 
obtained  by  RSS’ing  the  response  acceleration  magnitudes.  Table  2  summarizes  the  result  for 
each  SEA  subsystem.  For  example,  for  the  Top  Plate,  the  modal  transient  method  produced  3Ag 
and  -3.6g  temporal  peak  accelerations.  Using  PR=3,  maximum  global  response  magnitude  of 
9.0g  was  obtained.  PR=5  produced  5.4g  and  PR=7  3.9g  which  is  almost  identical  to  the  value 
obtained  from  the  modal  transient  method.  It  has  been  noted  that  PR=3  would  produce 
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conservative  results  enveloping  temporal  peak  accelerations  obtained  from  the  modal  transient 
method. 


Figure  5:  Response  of  the  Top  Plate  measured  at  Node  12359  using  modal  transient  method. 


Frequency,  Hz 


Figure  6:  Responses  of  the  3  SEA  subsystems  recovered  in  the  SRS. 
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Frequency,  Hz 


Figure  7:  Estimated  response  acceleration  magnitudes  in  the  time  domain  (Top  Plate). 


Table  2:  Maximum  global  response  magnitudes  for  SEA  subsystems. 
Green  shows  the  perfect  match  with  the  modal  transient. 


Maximum  Side  Plate 
Global  Response 
Magnitude,  g 

Maximum  Top  Plate 
Global  Response 
Magnitude,  g 

Maximum  Beam 
Global  Response 
Magnitude,  g 


SEA  (PR=7} 

SEA  (PR=5} 

SEA  (PR=3) 

FEM 

15.7g 

22.Dg 

36.  Gg 

22.7g.-12.6g 

a.gg 

5.4g 

9.0g 

3.4g,  -3.6g 

1.5g 

2.1  g 

3.4g 

2.6g,  -Z.7g 

Using  the  TF  of  the  SEA  model  and  the  input  acceleration,  response  waveforms  have  been 
calculated  for  each  SEA  subsystem.  Table  3  shows  the  maximum  global  response  magnitudes 
obtained  from  the  response  waveforms.  The  maximum  global  response  magnitude  of  the 
response  waveform  is  simply  the  maximum  value  of  the  waveform,  therefore  it  can  be  directly 
compared  with  the  maximum  global  response  magnitude  obtained  from  the  SRS  energy 
constraint  input  condition.  The  maximum  global  response  magnitudes  of  the  response  waveforms 
were  lower  than  those  from  the  SRS  energy  constraint  condition  in  Table  2.  For  example,  for  the 
Top  Plate,  the  SRS  constraint  condition  produced  the  result  (3.9g)  matching  the  result  from  the 
modal  transient  method  with  PR=7.  However,  the  response  waveform  produced  the  maximum 
global  response  magnitudes  of  2.8g  and  -4.7g.  If  necessary,  they  could  be  further  refined  using 
the  Magnitude  Correction  Factor,  C=1.39.  The  response  waveforms  were  further  refined  using 
the  phase  perturbations  based  on  bending  wave  numbers  of  the  Small  Side  Plate,  Side  Plate,  and 
Top  Plate  subsystems,  and  longitudinal  wave  numbers  of  the  Beam  subsystem.  It  seems  that  the 
phase  perturbation  may  have  improved  the  result  for  the  Top  Plate  subsystem,  however  no 
meaningful  improvements  were  found  from  the  Side  Plate  and  the  Beam  subsystems  from  the 
maximum  global  response  magnitude  perspective. 
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Table  3:  Maximum  global  response  magnitudes  using  response  waveforms  and  phase  perturbations. 


SEA  SRS 
Constraint 

SEA  (Response 
Waveform} 

SEA  (Phase 
Perturbed} 

FEM 

Magnitude 
Correction 
Factor,  C 

Nfaximum  Side  Plate 
Global  Response 
Magnitude,  g 

22.  Og 

13.0g,  -17.2g 

13.4g,  -19.0g 

22.7g,  -12.6g 

1.69 

Maximum  Top  Plate 
Global  Response 
Magnitude,  g 

3.9g 

2.8g,  -t.7g 

3.6g,  -3.5g 

3.4g,  -3.©g 

1.39 

Maximum  Beam 
Global  Response 
Magnitude,  g 

2.1  g 

1.4g,-1.&g 

1.4g.  -1.6g 

2.6g,  -2.7g 

1.50 

Figure  8  shows  the  response  waveform  of  the  Side  Plate  subsystem  compared  to  the  modal 
transient  result  recovered  at  Node  38801.  The  maximum  global  response  magnitude  was 
corrected  with  C=1.69  and  the  phase  perturbation  was  applied  without  an  energy  weighting.  It 
appears  that  the  response  waveform  follows  the  modal  transient  result  reasonably  well,  except  at 
0.05  seconds,  the  SEA  result  appears  to  be  opposite  to  the  modal  transient  result.  This  may 
simply  be  related  to  the  formulation  of  the  phase  perturbation  methods.  Figure  9  shows  the 
response  waveform  of  the  Beam  subsystem  compared  to  the  modal  transient  result  recovered  at 
Node  1.  The  maximum  global  response  magnitude  was  corrected  with  C=1.50  and  no  phase 
perturbation  was  applied.  The  SEA  result  showed  a  good  agreement  with  the  modal  transient 
result,  except  the  modal  transient  method  was  able  to  produce  the  “ringing”  of  the  beam  at  0.05 
seconds.  The  proposed  method  is  based  on  producing  maximum  responses  using  the  SRS  as  the 
input  condition  and  the  response  waveforms  are  only  “attached”  to  the  maximum  global  response 
magnitudes,  the  response  waveforms  would  not  be  able  to  simulate  local  modes  such  as  this 
“ringing”  mode.  However,  in  general,  the  proposed  method  was  able  to  produce  a  time  signature 
matching  reasonably  well  with  the  modal  transient  result  for  each  subsystem. 


is 


o  Phase  Perturbed  (SEA),  g 
o  NodeSSSOI  (FEM>,  g 


Time,  Seconds 


Figure  8:  Response  waveform  of  the  Side  Plate  subsystem  with  C=1.69  and  phase  perturbation. 
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Figure  9:  Response  waveform  of  the  Beam  subsystem  with  C=1.50  and  no  phase  perturbation. 

COMPLEX  VALIDATION  MODEL 

In  order  to  validate  the  proposed  pyrotechnic  shock  prediction  method,  a  complex  validation 
model  was  developed.  Figure  10  shows  the  complex  validation  model  undergoing  pyrotechnic 
shock  tests  using  the  Universal  Pyrotechnic  Shock  Simulator^^.  Figure  10  also  shows  the  SEA 
model  of  the  complex  validation  model.  In  this  paper,  only  the  preliminary  results  of  the 
correlation  study  are  provided.  More  than  50  pyrotechnic  shock  tests  have  been  conducted  to 
produce  transmissibility  functions  between  the  input  and  the  responses.  The  model  is  made  of 
aluminum  and  consists  of  plate  and  beam  like  elements  with  complex  geometries  and  various 
junctions  including  welding  and  multiple  fastener  combinations.  During  pyrotechnic  shock  tests, 
the  model  was  instrumented  with  more  than  15  Dytran  shock  accelerometers  to  measure 
response  time  histories  from  each  structural  member.  Dataflex-IOOOA  by  DSPCon  was  used  to 
measure  the  time  histories  with  100,000  Hz  sample  rate.  Figure  11  shows  the  SRS  input 
developed  from  the  temporal  accelerations  measured  at  the  system  interface  (Chi  and  Ch2).  The 
SRS  input  was  applied  as  the  energy  constraint  condition  at  the  interface  subsystem  of  the  SEA 
model.  Figures  12  and  13  show  the  response  spectra  recovered  at  Ch5  and  Chl3,  respectively.  It 
appears  that  the  SEA  results  are  well  correlated  with  the  test  data  and  the  predicted  SRS  levels 
from  the  SEA  model  successfully  envelope  the  test  data. 
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Figure  10:  Complex  validation  model  fully  instrumented  for  pyrotechnic  shock  tests  (L)  and  the  SEA 
model  (R).  More  than  50  pyrotechnic  shock  tests  were  conducted  to  develop  transmissibility  functions 
between  the  input  and  the  responses.  Maxi-Max  SRS  has  been  calculated  from  Chi  and  Ch2.  Green 
hourglass  represents  the  SRS  energy  constraint  input  condition  to  the  interface  SEA  subsystem. 
SEA-h  by  InterAC/LMS  was  used  for  the  SEA  modeling. 


Figure  11:  Maxi-Max  SRS  input  (Q=10)  calculated  from  Chi  and  Ch2. 
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Figure  12:  Comparison  of  the  response  spectra  recovered  from  the  proposed  SEA  method 

and  from  the  test  data  (Ch5). 


Figure  13:  Comparison  of  the  response  spectra  recovered  from  the  proposed  SEA  method 

and  from  the  test  data  (Chi 3). 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  proposed  pyrotechnic  shock  analysis  method  using  the  traditional  SEA  modeling  produced 
an  excellent  correlation  with  the  modal  transient  analysis  for  the  numerical  validation  model. 
Using  the  SRS  as  an  energy  constraint  condition,  the  proposed  method  was  able  to  estimate  the 
maximum  global  response  magnitudes  in  the  time  domain  and  produced  response  waveforms 
well  correlated  with  the  modal  transient  results.  In  addition,  preliminary  study  indicated  a  good 
correlation  with  the  test  data  measured  from  the  complex  validation  model.  Currently,  the  test 
data  from  the  complex  validation  model  are  being  investigated  for  each  structural  member  to 
further  validate  the  proposed  method.  The  proposed  method  may  be  used  to  produce  pyrotechnic 
shock  design  and  test  requirements  for  components  and  subsystems  due  to  the  SRS  requirement 
at  the  system  interface.  The  results  could  be  further  refined  with  response  waveform  predictions 
based  on  inverse  transform  of  the  product  of  the  input  acceleration  and  the  TF  of  the  SEA  model. 
The  response  waveform  may  be  used  to  estimate  various  engineering  variables  such  as  stress  and 
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velocity,  and  could  be  used  as  an  input  acceleration  to  shock  susceptible  components  for  the 
component  level  modal  transient  analysis. 
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